(* The following is to demonstrate the Monte Carlo method. Do not use for approximating value of the integral. *)
ClearAll[counter,n,x,y,f];
f[x_]:=E^(x^2)
NIntegrate[f[x],{x,0,1}]
counter[0]=0;n=5;
Table[{x=Random[Real,{0,1}],y=Random[Real,{0,E}],If[y≤f[x],counter[i]=counter[i-1]+1,counter[i]=counter[i-1]],
Print[" x=",x," y=", y," f(x)=", f[x], "  count=", counter[i], " approx area=" ,N[(counter[i]/n)*E]]},{i,1,n}];

(* The following is to demonstrate the Monte Carlo method. Do not use for approximating value of the integral. *)
ClearAll[counter,n,x,y,f];
f[x_]:=E^(x^2)
g1=Plot[f[x],{x,0,1},PlotRange->{{-0.2,1.2},{-0.5,3}},PlotStyle->Red];
NIntegrate[f[x],{x,0,1}]
counter[0]=0;n=5;
Table[{x[i]=Random[Real,{0,1}],y[i]=Random[Real,{0,E}],If[y[i]≤f[x[i]],counter[i]=counter[i-1]+1,counter[i]=counter[i-1]],
Print[" x=",x[i]," y=", y[i]," f(x)=", f[x[i]], "  count=", counter[i], " approx area=" ,N[(counter[i]/n)*E]]},{i,1,n}];
pts=Table[{x[i],y[i]},{i,1,n}];
g2=ListPlot[pts,PlotStyle->PointSize[0.015]];
g3=Graphics[Line[{{0,0},{1,0},{1,E},{0,E},{0,0}}]];
Show[g1,g2,g3]

(* The NIntegrate command is used to determine accuracy of the Monte Carlo method. Do not use it for approximating value of the integral using the Monte Carlo method. *)
ClearAll[counter,n,x,y,f];
f[x_]:=E^(x^2)
NIntegrate[f[x],{x,0,1}]
counter[n_][0]=0;n=10000;
Table[{x=Random[Real,{0,1}],y=Random[Real,{0,E}],If[y≤f[x],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n)*E ]]

(* The NIntegrate command is used to determine accuracy of the Monte Carlo method. Do not use it for approximating value of the integral using the Monte Carlo method. *)ClearAll[counter,n,x,y,f];
f[x_]:=E^(x^2)
g1=Plot[f[x],{x,0,1},PlotRange->{{-0.2,1.2},{-0.5,3}},PlotStyle->Red];
NIntegrate[f[x],{x,0,1}]
counter[n_][0]=0;n=10000;
Table[{x[i]=Random[Real,{0,1}],y[i]=Random[Real,{0,E}],If[y[i]≤f[x[i]],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n)*E ]]
pts=Table[{x[i],y[i]},{i,1,n}];
g2=ListPlot[pts,PlotStyle->PointSize[0.0001]];
g3=Graphics[Line[{{0,0},{1,0},{1,E},{0,E},{0,0}}]];
Show[g1,g2,g3]

(* The NIntegrate command is used to determine accuracy of the Monte Carlo method. Do not use it for approximating value of the integral using the Monte Carlo method. *)ClearAll[counter,n,x,y,f];
f[x_,y_]:=E^(-x^2)E^(-y^2)
Plot3D[f[x,y],{x,0,1},{y,0,1}]
NIntegrate[f[x,y],{x,0,1},{y,0,1}]
counter[n_][0]=0;n=1000;
Table[{x=Random[Real,{0,1}],y=Random[Real,{0,1}],z=Random[Real,{0,1}],If[z≤f[x,y],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n) ]]

(* The NIntegrate command is used to determine accuracy of the Monte Carlo method. Do not use it for approximating value of the integral using the Monte Carlo method. *)ClearAll[counter,n,x,y,z,f];
f[x_,y_]:=E^(-x^2)E^(-y^2)
g1=Plot3D[f[x,y],{x,0,1},{y,0,1}];
NIntegrate[f[x,y],{x,0,1},{y,0,1}]
counter[n_][0]=0;n=100;
Table[{x[i]=Random[Real,{0,1}],y[i]=Random[Real,{0,1}],z[i]=Random[Real,{0,1}],If[z[i]≤f[x[i],y[i]],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n) ]]
pts=Table[{x[i],y[i],z[i]},{i,1,n}];
g2=Graphics3D[Point[pts],PlotStyle->PointSize[0.0001]];
Show[g1,g2]

ClearAll[counter,n,x,y,z,f];
f[x_,y_]:=Sqrt[1-x^2-y^2]
Plot3D[f[x,y],{x,0,1},{y,0,1}]
NIntegrate[f[x,y],{x,0,1},{y,0,Sqrt[1-x^2]}]
counter[n_][0]=0;n=1000;
Table[{x=Random[Real,{0,1}],y=Random[Real,{0,1}],z=Random[Real,{0,1}],If[f[x,y]∈Reals && z≤f[x,y],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n) ]]

ClearAll[counter,n,x,y,z,f];
f[x_,y_]:=Sqrt[1-x^2-y^2]
g1=Plot3D[f[x,y],{x,0,1},{y,0,1}];
NIntegrate[f[x,y],{x,0,1},{y,0,Sqrt[1-x^2]}]
counter[n_][0]=0;n=100;
Table[{x[i]=Random[Real,{0,1}],y[i]=Random[Real,{0,1}],z[i]=Random[Real,{0,1}],If[f[x[i],y[i]]∈Reals && z[i]≤f[x[i],y[i]],counter[n][i]=counter[n][i-1]+1,counter[n][i]=counter[n][i-1]]},{i,1,n}];
Print[N[(counter[n][n]/n) ]]
pts=Table[{x[i],y[i],z[i]},{i,1,n}];
g2=Graphics3D[Point[pts],PlotStyle->PointSize[0.0001]];
Show[g1,g2]